Multifractality approach of a generalized Shannon index in financial time series

Multifractality is a concept that extends locally the usual ideas of fractality in a system. Nevertheless, the multifractal approaches used lack a multifractal dimension tied to an entropy index like the Shannon index. This paper introduces a generalized Shannon index (GSI) and demonstrates its application in understanding system fluctuations. To this end, traditional multifractality approaches are explained. Then, using the temporal Theil scaling and the diffusive trajectory algorithm, the GSI and its partition function are defined. Next, the multifractal exponent of the GSI is derived from the partition function, establishing a connection between the temporal Theil scaling exponent and the generalized Hurst exponent. Finally, this relationship is verified in a fractional Brownian motion and applied to financial time series. In fact, this leads us to proposing an approximation called local fractional Brownian motion approximation, where multifractal systems are viewed as a local superposition of distinct fractional Brownian motions with varying monofractal exponents. Also, we furnish an algorithm for identifying the optimal q-th moment of the probability distribution associated with an empirical time series to enhance the accuracy of generalized Hurst exponent estimation.


Introduction
The scaling property of a mathematical function is a crucial tool for understanding system variations based on exponents associated with universality classes [1].In stochastic processes, this property is explored by examining the invariance of the probability distribution pattern of a random variable [2].For stationary time series, the identification of a single exponent, denoted as δ and referred to as the scaling exponent, is found to be adequate.Early pioneers such as B. Mandelbrot and H. Stanley played significant roles in calculating this exponent in time series [1,[3][4][5][6].Since the 1980s, various methods, including the precise Shannon entropy calculation in generated sub-series via the diffusive trajectory algorithm, have been proposed for determining the scaling exponent δ [7].Notably, the diffusive trajectory algorithm generates multiple sub-series by adding consecutive terms, akin to a simple Brownian motion process, and computes the theoretical scaling exponent of a Le ´vy flight [7,8].However, challenges arise as the scaling exponent δ doesn't always coincide with the Hurst exponent (H) which describes the evolution of the second moment around the origin in an anomalous diffusive process [9][10][11][12].
On the other hand, the Theil index (T) is an inequality measure devised by economist Henri Theil and formulated in terms of an entropy index [13].Also, this index boasts a decomposition property wherein the overall income inequality of a population is derived from subgroup inequalities [14,15].Employed in econophysics, the Theil index has been used to explore correlations in time series, transforming them into entropy time series based on time window size [16].Its applications extend to comprehending equilibrium states in free market models [17], analyzing regional changes in foreign aid distribution using an entropic approach [18], and studying income distribution in countries [19].The Theil index also finds relevance in characterizing racial segregation, measuring redundancy in data, and assessing diversity [20][21][22].
While comparing Theil indices across populations of different sizes poses challenges, normalization by the logarithm of data number renders them comparable.Normalized, the Theil index mirrors the Gini index, introduced by sociologist Corrado Gini in 1912, with both reaching 0 for maximum equality and 1 for maximum inequality [23].Furthermore, for a parametric family of probability distributions, the Gini index is considered a measure of dispersion linked to the variance (X 2 ) [24].Similarly, the Theil index for major parametric income distributions is expressed in terms of variance, indicating a distribution-dependent relationship between these two quantities [25].Thus, the Theil index, often associated with variance [26], is seen as a measure of dispersion.Consequently, considering the empirical relationship between the variance X 2 and the mean (M 1 ), that is called temporal fluctuation scaling, an expected relationship between T and M 1 emerges, as established in the literature and termed temporal Theil scaling [27].In fact, it has also been shown that the temporal fluctuation scaling is not always satisfied in the diffusive trajectory algorithm [27].Therefore, in complement to the previous ideas, the inquiry persists regarding the possibility of establishing a relationship between the Hurst exponent and the temporal Theil scaling exponent within the context of the diffusive trajectory algorithm.
Regarding the temporal Theil scaling, it is important to mention that this temporal scaling is described in terms of the Shannon index and the diffusive trajectory algorithm [27].Thus, remembering that the Shannon index is described in terms of the Theil index [13], and this in turn in terms of the generalized entropy index GE(q) where q is associated to q-th moment of the probability distribution of income of a population [20], the question naturally arises as to how to extend the Shannon index so that it depends on the q-th moment of the probability distribution since this opens the possibility of exploring new types of temporal scaling.Nevertheless, this must be done carefully, keeping in mind that the generalized entropy index has a removable discontinuity for q = 0 and q = 1.Also, it is observed that in principle the only requirement to calculate this Shannon index is that the probability distribution has a support of positive real values as it happens for a positive definite measure as a probability measure.
Multifractality is a concept that extends traditional ideas of fractality locally within a system F [28], originated as a neologism coined by Frish and Parisi in their work on turbulence in 1983 [29,30].This term was further affirmed by Mandelbrot [30].Over time, multifractality has been explored through two prominent approaches: the Structure Functions Approach (SFA) and the Partition Function Approach (PFA).In the case of the SFA, this is proposed by Kantelhardt et al. in 2002 [31], defining the structure functions in terms of the q-th moment of the probability distribution associated with time-lagged increments δ in the system F.Then, a power law relationship is then assumed for these structure functions to link the exponents to the generalized Hurst exponent H(q), generalized fractal dimension D q , or the mass exponent function τ(q) [31,32].
In contrast, the PFA, introduced by Grassberger and Procaccia in 1983 [33], defines a partition function in terms of a probability measure f F within the system F, akin to the usual ensembles of statistical mechanics.The behavior of the system is characterized by calculating the generalized fractal dimension in terms of the Renyi entropy [33].In fact, q = 0, q = 1, and q = 2 correspond to the usual fractal dimension [34], the dimension of information calculated with Shannon's entropy [35], and the dimension of correlation [33], respectively.
Furthermore, several studies in the literature have concentrated on examining and elucidating certain properties regarding the potential relationship between entropy measures and long-range memory measures, either implicitly through the auto-correlation function [36][37][38][39] or explicitly via the Hurst exponent [40].Notably, the probability of a particular pattern p π (d), characterized by a delay d and a permutation π, being reiterated within a time series enables the establishment of properties akin to the temporal correlation of the series in a much-simplified way [36].For instance, analytical expressions are derived for fractional Brownian motion and processes with stationary Gaussian increments, where the auto-correlation function hinges solely upon the probabilities of patterns of order p π (d = 2) [36,37].Nevertheless, pattern order analysis encounters complexities, particularly when the pattern's order surpasses four, leading to intricate or dependent on many parameters at the computational level functions [36,37].
Moreover, the generalized fractal dimension of a system serves as a direct measure of the system's entropy [38,40], suggesting an inherent connection between Renyi entropy and the generalized fractal dimension [33,39,40].The above has motivated the interest in employing entropic analyses, such as multi-scale entropy analysis, to establish a direct linkage between the Hurst exponent of fractional Brownian motion and the sample entropy (S E ), assuming that S E conforms to a q-exponential function [40].Yet, to obtain a broader perspective of the relationship between entropy as an indirect measure of a system's fractality, it is imperative to recognize several limitations in the existing literature.
Firstly, the most common entropy measures, like Renyi entropy, are non-extensive quantities in the physical realm.This non-extensivity implies that the total entropy of a system is not computed as the sum of individual entropies when the system is divided into several regions.Consequently, entropy is not a scale-invariant or self-similar measure, which is fundamental in a fractal system.Secondly, while the diffusive trajectory algorithm measures the Hurst exponent of a fractional Brownian motion exactly [7], it prompts questions regarding the feasibility of extending such an analogy to systems that do not satisfy a fractional Brownian motion.Thus, the possibilities of extending these approaches to normalized measures, such as the Shannon index, are raised.Also, it is highlighted that this extension would arise naturally with the temporal Theil scaling exponent as the latter is linked to measures such as the Shannon index and the diffusive trajectory algorithm [27].Accordingly, to establish a clearer relationship of entropy as a measure of the fractality of a system, one solution is to establish a connection between the generalized Hurst exponent and the temporal Theil scaling exponent as addressed in this work.
Therefore, in this paper, a review of the theoretical approaches used to address multifractality is made.From the above, an extension of the Shannon index is proposed, and then a novel theoretical relationship is deduced between the generalized Hurst exponent with a newly introduced multifractal exponent termed the multifractal exponent of the generalized Shannon index.To do this, the concept of multifractality within a system F is defined in Section 2, delving into two schemes for computing multifractal exponents: the structure function approach and the partition function approach.Actually, one notable multifractal exponent is the generalized Hurst exponent H(q), where q represents the q-th moment of the associated probability distribution with a measurement on the system F. Later, in section 4 the extension of the Shannon index taking into account the generalized entropy index is made.Thus, this extension is given the name generalized Shannon index and allows to define a newly multifractal exponent β TTS (q) in such a way that the generalized Hurst exponent is related to the temporal Theil scaling exponent as shown in Eq (31) from the subsection 5.In Section 6, the proposed relation is validated for fractional Brownian motion, and Section 7 demonstrates its application to a two financial time series, offering a optimal selection algorithm for the most significant q-th moment in estimating the generalized Hurst exponent.

Multifractality in complex systems
Fractality typically pertains to the geometric attributes of an object, emphasizing properties like self-similarity, scale invariance, and nowhere differentiable, with its fractal dimension D surpassing the topological dimension [5,41].In other words, a fractal object exhibits recurring patterns across scales, featuring irregularities that lack smoothness in terms of differential calculus.Also, its length, area, or volume scales following a power law, where the exponent is a non-integer value surpassing the conventional dimension [5,41].Indeed, the characterization of the fractal dimension D typically employs the box counting method.This approach involves counting the units necessary to cover the fractal set denoted by N, considering a unit size δ > 0 as the measuring scale [34].Thus, assuming [34] N � d À D ; ð1Þ it follows that the fractal dimension is the slope of a log-log graph of the number of units needed to cover F as a function of the unit of measure δ.It is important to mention that there are other measures to describe fractality and that they have emerged in different branches of science such as the information dimension [35], the correlation dimension [33], the Lyapunov dimension [42], the Higuchi dimension [43], or the Haussdorf dimension [44].
Although, there are fractal systems where a single fractal dimension is not enough to describe the collective behavior of the whole system and in such a case the concept of multifractality is reached.Specifically, a system F is considered multifractal if given a measure f F on the system it is satisfied that where x is a specific instant of time t or spatial pointr under observation, ||δ|| > 0, and y(x) is called the singularity exponent [28].Thence, multifractality is a property wherein a system exhibits local fractal behavior based on the specific instant of time t or spatial pointr under observation [28].Historically, the term multifractality was first coined by Frish and Parisi in their 1983 work on turbulence [29,30], a concept later confirmed by Mandelbrot [30].Concurrently, in the same year, Grassberger and Procaccia introduced the concept, utilizing Renyi entropy [33].Presently, two predominant approaches, known as the partition function approach (PFA) and the structure functions approach (SFA), are employed to explore and understand multifractality.
In the SFA, proposed by Kantelhardt et.al in 2002 [31], it is found that the increments of the measure f F characterize the entire behavior of the system, i.e, estimating the q-th moment of the probability distribution associated with the measure f F is satisfied for all δ > 0 that where H(q) is the generalized Hurst exponent [31,32].Furthermore, when H(q) is constant, the system is said to have monofractal behavior and H(q) = H corresponds to the usual Hurst exponent, whereas if qH(q) is a nonlinear function then the system has multifractal behavior and is a strong argument against the Brownian, Fractional Brownian, Le ´vy and Fractional Le ´vy models, which are additive models, so qH(q) are portions of straight lines [32].It goes without saying that it is also a known fact that a possible explanation for the origin of multifractality is using geometric Tweedie models [45,46].On the other hand, in the PFA using the same idea of box counting in fractality, it is assumed that the measure f F is a probability measure such that if the system F is divided into boxes of size δ > 0 the following two equalities are satisfied where Bðd; sÞ is the s-th box of δ > 0 size in the system.Thus, the partition function O F ðq; dÞ is defined by [28,33] O F ðq; dÞ ¼ Thence, if the Renyi entropy Iðq; dÞ is defined as [33] Iðq; dÞ ¼ lim and the generalized fractal dimension D q is [28, 33] Iðq; dÞ An important fact to note here is that the generalized fractal dimension with q = 0, q = 1, and q = 2 corresponds to the usual fractal dimension [34], the information dimension calculated with Shannon entropy [35], and the correlation dimension [33], respectively.Therefore, from the generalized fractal dimension D q , the generalized Hurst exponent H(q) is defined as [28,33] Finally, to complete the whole panorama of multifractality, it is crucial to highlight two additional functions that gauge the degree of multifractality within the system.These are referred to as the mass exponent function τ(q) and the singularity spectrum ϑ(y).The mass exponent function is related with the generalized Hurst exponent as follows [28,33] tðqÞ ¼ qHðqÞ À 1: ð10Þ In the same way, the singularity spectrum is expressed in terms of the mass exponent function through the Legendre transform, that is, Indeed, the singularity spectrum ϑ(y) corresponds to a concave function that gathers all the singularity exponents or singularity strengths y(t) of the system [28,33].

Multifractal Detrended Fluctuation Analysis
Now, to grasp the computation of the generalized Hurst exponent in empirical time series data, the Multifractal Detrended Fluctuation Analysis (MF-DFA) method is presented [28,31,32].It is essential to note that various methods exist for calculating the generalized Hurst exponent, including the detrended moving average (DMA) [47], geometric method-based procedures (GM algorithms) [9], absolute value method (AVE) [10], fractal dimension algorithms (FD) [48], Generalized Hurst Exponent (GHE) [49], triangle total areas (TTA) [12], and the KS method [50].Nevertheless, preference is given over the MF-DFA considering that for this generalized Hurst exponent method its canonical measure function is known [51].
Multifractal Detrended Fluctuation Analysis (MF-DFA) is a method designed on innovation time series such that the generalized Hurst exponent H(q) is computed over a time series fZðtÞg N t¼1 with the following steps [31]: 1. Determine the profile of the time series as 2. Divide the profile U(t) into N δ = bN/δc non-overlapping segments of length δ > 0. Thus, if the time series is not a multiple of the considered time scale δ > 0, then the same procedure is repeated over time series starting from the opposite end.
3. Calculate the local trend Ũ for each of the N δ segments (2N δ if N=d = 2 N) using some approximation technique, such as a detrended fluctuation analysis (DFA or polynomial fitting), or a detrended moving average (DMA).The detrended residuals f�ðtÞg N t¼1 are defined by 4. Estimate the variance for each of the N δ segments (2N δ if N=d = 2 N) as for the s-th segment of size δ > 0. It is important to note that in the second line of Eq (14), there is the restriction s > N δ .
5. Average over all segments to obtain the q-th order overall detrended fluctuation function 6. Repeat steps 2 to 5 for several time scales δ.
7. Determine the scaling behavior of the q-th order overall detrended fluctuation by analyzing log-log plots Fðq; dÞ versus δ > 0 for each value of q 2 R. If the time series fZðtÞg N t¼1 are long-range power-law correlated, Fðq; dÞ increases, for large values of δ > 0, as a power-law [31], Fðq; dÞ � d HðqÞ ; ð16Þ where H(q) is the generalized Hurst exponent.
Here it is worth mentioning that for values q > 0, the behavior of the generalized Hurst exponent is dominated by the segments with the highest variance F2 (δ, s) (see Eq ( 14)), which implies that H(q) describes the scaling behavior of segments with large fluctuations.In the same way, for values q < 0, the generalized Hurst exponent H(q) explains the scaling behavior of segments with small fluctuations, which are usually characterized by a higher exponent [28,31].Now, if the form of the expressions ( 6) and ( 15) is compared, a high similarity is observed to the point that it is affirmed that the partition function of the MF-DFA method is Indeed, defining a canonical measure μ(q, δ, s) = F q (δ, s)/O MF (q, δ) in terms of this partition function and a simple moving average, the generalized fractal dimension D q (see Eq (8)), the singularity strength yðqÞ ¼ dt dq , and the singularity spectrum ϑ(y) (see Eq (11)) are calculated [51].Also, it is worth noting that many times the fluctuations or variances that are observed locally in the residuals of the system, measured through F 2 (δ, s), do not exceed the value of 1, which implies that the partition function (17) satisfies the expressions (4) and ( 5) as required for the PFA.
Finally, this section is concluded by presenting a stochastic process that is defined with the classical Hurst exponent called fractional brownian motion (fBm).Fractional Brownian motion is the generalization of Brownian motion and first introduced in 1953 by P. Le ´vy although the name was actually given by B. Mandelbrot once he recognized its importance [52,53].Actually, the fractional Brownian motion with Hurst exponent H = H(q = 2), denoted by {B H (t)} t�0 , is characterized by the following four properties [52]: 4. It has a Gaussian distribution for its stationary increments, that is, , for all t � s � 0, where N ðm; sÞ represents a normal distribution with mean μ and standard deviation σ.
It is noteworthy that, in a historical context, fractional Brownian motion was utilized even before P. Le ´vy's contributions. A. N. Kolmogorov employed fractional Brownian motion in 1940 to investigate spiral curves in Hilbert space [54].Furthermore, its applications extend to the realm of the random Fourier transform [55], and the exploration of correlations in stochastic processes with stationary increments of order n [56].The significance of fractional Brownian motion grew over time, culminating in its formal definition by B. Mandelbrot and J. Van Ness in 1968, encapsulated in the stochastic fractional equation [52,53] where ξ(t) is a white noise and t a I n t b is the right-handed-sided Riemann-Liouville fractional integral of order H þ 1 2 > 0. Thus, if H = 1/2, the fractional Brownian motion is an usual Brownian motion, if H > 1/2, the increments of the process are positively correlated, and if H < 1/2, the increments of the process are negatively correlated.Here it is worth mentioning that the Hurst index is a measure of long-range memory in time series.Thus, it is related to the auto-correlation of the time series and the rate at which it decreases as the delay between pairs of values increases.

Generalized Shannon index
In the preceding section, multifractality was identified as an expansion of the fractal characteristics of a system F, wherein the fractal behavior is contingent upon the observation point, specifically on the neighborhood in which it evaluates the singularity exponent y(t).Thus, various multifractality concepts involve the definition of limits and, consequently, of neighborhoods around a point (refer to Eqs (7), ( 8) and ( 9)).Thence, to delve into the theoretical understanding of the relationship with the temporal Theil scaling exponent α TTS (t), it becomes imperative to extend the Shannon index through the utilization of the generalized entropy index GE(λ).This extension broadens the Theil index's definition and the foundations of temporal Theil scaling.Now, we define the generalized entropy index GE(λ) as [20] where N represents the size of the population, y i is the income for case i, and � y is the sample average.Hence, with small values of λ, the generalized entropy index GE(λ) exhibits sensitivity to minor incomes, while larger λ values render the index sensitive to major incomes.Notably, when λ = 1, the generalized entropy index GE(λ) corresponds to the Theil index T. Also, the Shannon index T was defined as a normalization of the Theil index [57].Therefore, if it is observed that this normalization value is made with respect to the maximum value that the Theil index can take, then the generalized Shannon index SðlÞ is defined as Thus, changing the dummy variable λ to q and defining the normalized income x k ¼ y k = P N j¼1 y j , it is clear that the generalized Shannon index is Note that Eq (21), due to the way it was defined, is maximum when there is maximum equality in the income distribution, that is, x k = 1/N, for all k 2 {1, 2, . .., N}.From this, taking a time series of positive values fZ t g N t¼1 , the diffusive trajectory algorithm generates multiple time sub-series X(t, s), such that the generalized Shannon index at time t is where Xðt; sÞ is the normalized s-th diffusive trajectory defined as In Eq (23), Xðt; sÞ ¼ P t i¼1 Z iþs is known as the s-th diffusive trajectory of fZ t g N t¼1 .Finally, the temporal Theil scaling is defined as the following power law between the normalized Shannon index and the normalized mean of each diffusive trajectory [27] Sð1; tÞ where M 1 (t), K TTS (t) and α TTS (t) are the mean of the s-th diffusive trajectory, the constant and exponent of the temporal Theil scaling, respectively.Additionally, 5 Temporal Theil scaling exponent as multifractal exponent Now, with the extension made to the Shannon index T , notice that Eq (21) is similar to Renyi entropy (see Eq (7)) in the sense that both are measures of entropy over an arbitrary complex system or data set as a time series.Besides, the Renyi entropy is not an extensive quantity in the physical sense, which implies that paradoxes such as the Gibbs paradox are generated [58], while the generalized Shannon index is a homogeneous function of degree q, that is, if all the values in the system increase by an amount β > 0 it follows that since the maximum value is now proportional to β q−1 N q−1 , while the normalized income does not change under scaling.Thus, the generalized Shannon index is a more precise quantity to use in describing the scale invariance of a system as well as not having the Gibbs paradox problem when q = 1.Therefore, it is useful to be able to extend the ideas of multifractality by means of the generalized Shannon index in such a way that if the expressions ( 6), ( 7), ( 17) and ( 22) are compared, the partition function of the generalized Shannon index is defined by analogy as It is important to note that the partition function ( 27) has a functional form already similar to that of the temporal Theil scaling.Although, a last step is missing regarding normalizing the generalized Shannon index with respect to the maximum value in time, that is, S max ðqÞ.
Returning to the Eqs ( 22) and ( 27), it is clear that any type of normalization with respect to the maximum in the partition function of the generalized Shannon index depends on the size of the system, that is, on the length of the time series, and in this way it is assumed in forward that Sðq; dÞ where δ max denotes the time step where generalized Shannon index is maximum, M q (δ) is the q-th moment of the diffusive trajectories time series at time δ, M q,max = max δ {M q (δ)}, D 0 is the fractal dimension of the system or usual dimension, λ TTS−MF > 0 is the coupling parameter of the generalized Hurst exponent with the temporal Theil scaling exponent, and k TTSÀ MF 2 R is the coupling parameter of the usual dimension with the temporal Theil scaling exponent.Also, note that these coupling parameters arise as an effect of the size of the system, that is, the length of the time series N. Therefore, it is now possible to define a new multifractal exponent based on generalized Shannon index partition function and in such a way that the temporal Theil scaling exponent corresponds to a particular case of it.Thus, the multifractal exponent of the generalized Shannon index is defined as where is the critical parameter associated with an ordered phase in a condensed matter phase transition [27,59,60].Observe that the limit of Eq (29) is with respect to time δ !δ max which is well defined since the time instant in which the maximum occurs for the generalized Shannon index is the same for the moments around the origin of the distribution.In fact, said instant of time, as seen at the beginning of this section, corresponds to the moment in which the diffusive trajectory has maximum equality in analogy to an income distribution, which implies the lowest generalized entropy index over time for the multiple diffusive paths.Finally, if Eq (29) is applied for q = 1, the temporal Theil scaling exponent in terms of the generalized Hurst exponent H(q) is where the last step is only valid in monofractal systems such as the fractional Brownian movement with Hurst exponent H.

Temporal Theil scaling in fractional Brownian motion
This section verifies the validity of the generalized Shannon index multifractality approach proposed in the previous section.To do this, starting from Eq (31), the simulation of multiple fractional Brownian motion trajectories is performed and its exponent of the temporal Theil scaling is calculated to corroborate that the found expression works in this type of stochastic process.With the above, the process to follow to corroborate Eq (31) is: 1. Simulate N fBm fractional Brownian motion trajectories by setting T fBm time steps for different values of the Hurst exponent H.For this, the following aspects are taken into account: • The distribution of increments of a fractional Brownian motion, i.e. g H (t) = B H (t + 1) − B H (t), is called fractional Gaussian noise, and the latter has already been computationally implemented quite efficiently in Python [61].Indeed, the implementation process outlined in Python by Rydin Gorjao et al. follows a rough procedure [61].Initially, a circulant covariance matrix is defined, with its entries determined by the covariance function given by gðkÞ , for all k � 0. Subsequently, this matrix is diagonalized.Finally, two identically distributed standard normal random variables are generated, and the fast Fourier transform is applied in conjunction with the eigenvalues and eigenvectors of the circulant covariance matrix [62].This process yields a sample of fractional Brownian motion by utilizing the accumulated increments of the fractional Gaussian noise sample mentioned earlier.Notably, this technique is recognized in the literature as the Wood-Chan or Davies-Harte method [62][63][64].
• The estimate of the Hurst exponent in short time series has a high uncertainty as it has already been shown in [65].Then, a selection criterion is generated on each fractional Brownian motion path such that if the computationally found Hurst exponent, say H exp , is in the range [H − ε fBm , H + ε fBm ], then the simulation is a valid simulation.Otherwise, the process is repeated until completing the N fBm simulations.Here it is worth mentioning that once the sample of a fractional Brownian motion with Hurst exponent H is generated in the previous step, the computational Hurst exponent H exp is estimated in each simulation through the MF-DFA method mentioned in the section 3 and implemented in Python by Rydin Gorjao et.al. through a function called MFDFA [61].
2. Estimate the increments of the fractional Brownian motion sample (fBm), the absolute value of the fractional Brownian motion, and the volatilities of the fractional Brownian motion for each simulation using the following definitions ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi respectively.It is important to mention that in Eq (34), σ[�] represents the standard deviation of the random variable.Furthermore, the term inside the square root is usually known as standard score or z-score.Note that the expected value E½�� or the standard deviation σ[�] used in the definition of the z-score are statistical properties that are obtained at each time step from the usual formulas for the sample mean and the unbiased standard deviation of a sample, respectively.Thence, the total time series is tripled.
3. Apply the diffusive algorithm described by Xðt; sÞ ¼ P t i¼1 B H ði þ sÞ, for the original time series of fractional Brownian motion, absolute value of the fractional Brownian motion, and volatilities of the fractional Brownian motion.

Calculate the cumulative mean, and the Shannon index in each of the diffusive trajectory
time series by adding new data for each time t.(24) to find the temporal Theil scaling coefficient K TTS (t) and the temporal Theil scaling exponent α TTS (t).It is important to mention that in all the iterations of this process, one cost function is estimated on the adjustment of the temporal Theil scaling given by the coefficient of determination R 2 TTS .6. Filter the simulations where R 2 TTS exceeds a certain threshold value r fBm � 0 and average the values obtained for the temporal Theil scaling exponent in the different simulations of a fractional Brownian motion with Hurst exponent H. Thus, when repeating the process with different values of the Hurst exponent, we have the average of the temporal Theil scaling exponent α TTS (H) as a function of the Hurst exponent H. (31) to find the coupling parameter of the generalized Hurst exponent with the temporal Theil scaling exponent λ TTS−MF > 0 and the coupling parameter of the usual dimension with the temporal Theil scaling exponent k TTSÀ MF 2 R. It is important to mention that in all the iterations of this process, two cost function is estimated on the adjustment of Eq (31) given by the coefficient of determination R 2 and the p-th mean absolute error MAE p (n H ) defined by

Perform linear regression based on Eq
where n H is the number of the different values used for the Hurst exponent, and e k are the residuals obtained after adjusting for least squares.
Finally, before going on to the results after having done this process, it is worth mentioning that the parameters selected for the simulations were N fBm = 200, T fBm = 512, Δt = 0.001, ε fBm = 0.02, r fBm = 0.98, H 2 {0.4,0.46, 0.52, 0.58, 0.64, 0.7, 0.76, 0.82, 0.88} (then n H = 9).Regarding the values chosen for the Hurst exponent, it goes without saying that values are taken where the Hurst exponent is persistent (H ≳ 0.5) since the MF-DFA method for multifractality becomes imprecise for signals strongly anti-correlated when H(q) is close to zero.Also, taking into account some methods developed for the estimation of the Hurst exponent with short time series that use the entropy of Shannon, the values of H * 0.9 are limited to avoid the problem of bias of the diffusive trajectory algorithm [66].In addition, all the codes made are published in the Github [67].At last, to obtain the uncertainty in the parameters estimated by the least squares regression in some adjustments, the covariance matrix Cov(�) is calculated with the vector of residuals of the fitẽ, that is, Covð�Þ / ẽẽ T .Then, the square root of the respective term on the diagonal of the covariance matrix is taken as the error in the fit parameters.
Fig 1 illustrates the temporal Theil scaling exponent as a function of the Hurst exponent under the specified parameters.Notably, some datasets display a considerable standard deviation even after averaging multiple simulations, suggesting the potential for improvement through an increased number of simulations N fBm or by utilizing a fractional Brownian movement with a greater number of time steps T fBm .However, it is essential to acknowledge the computational challenges associated with the diffusive trajectory algorithm, whose execution time and memory usage grows proportionally to OðT 2 fBm Þ.In the regression analysis outlined in (31), the coefficient of determination R 2 for the time series of absolute value of the fractional Brownian motion and the time series of volatilities of the fractional Brownian motion are reported as 54.21% and 84.96%, respectively.This data, along with the coupling parameters obtained in each case (considering D 0 = 1 for time series), is presented in Table 1.In fact, it is observed that the coupling constant between the temporal Theil scaling exponent and the fractal dimension κ TTS−MF is negative in both cases, indicating an inversely proportional relationship between these exponents.Conversely, the coupling constant between the temporal Theil scaling exponent and the generalized Hurst exponent λ TTS−MF is positive but two orders of magnitude smaller, suggesting that the multifractal exponent of the generalized Shannon index is more sensitive to changes in the fractal dimension of the fractional Brownian motion.Finally, it is crucial to highlight that the temporal Theil scaling of the original time series is not performed, as a fractional Brownian motion is not positively defined by its positive or negative increments.
Therefore, the multifractal exponent of the generalized Shannon index demonstrates a high degree of precision for a fractional Brownian movement, as evidenced by the satisfaction of the expression (31), especially in the case of volatilities of the fractional Brownian motion, where the coefficient of determination is R 2 = 84.96%.However, the time series of absolute value of the fractional Brownian motion may appear to deviate from this trend, as indicated by the R 2 value of 54.21%, but in such a case it should be remembered that the small values close to zero associated with absolute value of the fractional Brownian motion magnify the impact of small fluctuations in the time series.To address this, it is recommended to compute the generalized Hurst exponent H(q) with q < 0, as opposed to the conventional Hurst exponent H (q = 2) = H.Furthermore, remembering that the volatilities of the fractional Brownian motion definition precisely mitigates the issue of small fluctuations when the z-score is calculated (see Eq (34)).

Application of the multifractal exponent of the generalized Shannon index
Finally, this section presents the association between the temporal Theil scaling exponent and the generalized Hurst exponent in empirical time series, a relationship previously substantiated in Section 4 and validated for fractional Brownian motion in Section 6. Thence, when dealing with time series of empirical data with an unknown Hurst exponent, a comprehensive exploration is required to obtain multiple temporal Theil scaling exponents and generalized Hurst exponents from the same time series.The proposed methodology involves accumulating data from the time series to calculate multiple temporal Theil exponents and generalized Hurst exponents, thereby constructing a substantial sample of potential values for these parameters.Thus, comparing α TTS and H(q) at the same instant of time, it is possible to observe the relationship between these two quantities.
For this purpose, the Taylor series of the multifractal exponent of the generalized Shannon index for a fractional Brownian motion with Hurst exponent H fBm around q = 1 is considered, Table 1.Parameters of the relationship between the multifractal exponent of the generalized Shannon index β TTS (q) and the generalized Hurst exponent H(q) together with its coefficient of determination R 2 obtained after the fit with the expression (31) for a fractional Brownian motion (fBm).

Time series
Coupling type Coupling constant value MAE 1 (×10 that is where is calculated by Eq (29), and the Taylor series is truncated to first order by the linear relationship that was verified for the fractional Brownian motion in the section 6.Now, equating (36) with (31), and solving for q − 1, we obtain The next step is to note that the generalized Hurst exponent H(q) in terms of the multifractal exponent of the generalized Shannon index β(q) satisfies that where the expression (29) has been used such that the multifractal exponent of the generalized Shannon index β(q) is proportional to the generalized Hurst exponent H(q).Thus, the n-th derivative of the generalized Hurst exponent, with n 2 N, is Then, the Taylor series of the generalized Hurst exponent around q = 1 is where δ j,0 is the Kronecker delta, β(1) = α TTS and the expressions (38) and (39) were used.
Furthermore, we define Consequently, motivated by Eq (37), it is assumed that which will be called hereinafter local approximation of fractional Brownian motion (LA-fBm).Thus, it is interpreted that in a small enough neighborhood of size ε > 0 around q, the generalized Hurst exponent behaves as a fractional Brownian motion with Hurst exponent HfBm ¼ H fBm ðq; εÞ.In other words, the generalized Hurst exponent is approximated by chunks of linear functions with different Hurst exponents, or equivalently, by multiple fractional Brownian motions.Consequently, the LA-fBm is an approximation in which a multifractal system is locally composed of different fractional Brownian motions such that some of them contribute more than others depending on the type of system and the time or space in which it is observed.It is important to note then that the factor that is preserved independent of LA-fBm is a TTS l À 1 TTSÀ MF , which implies that the expression ( 40) is written as where N corresponds to truncation up to the N-th order and R N (α TTS ) is a residual term or approximation error.Also, note that to guarantee convergence of Eq (43) it suffices that |α TTS − λ TTS−MF (H fBm (q, ε) + 1) + κ TTS−MF D 0 | < 1, as happens for a fractional Brownian motion (see Eq (31) with H = H fBm (q, �)), and hence the name of the approximation made.
Therefore, for simplicity, henceforth assume that where W is the degree of a polynomial of fit and fb n ðqÞ ¼ c n ðqÞl À n TTSÀ MF g W n¼1 the associated coefficients.Note that if W = 1, then it is possible to rewrite the expression (31) for a fractional Brownian motion as a linear function where b 1 ðqÞ ¼ l À 1 TTSÀ MF and b 0 ðqÞ ¼ D 0 k TTSÀ MF l À 1 TTSÀ MF À 1.Thus, the Eq (44) corresponds to considering the effects of multifractal time series as a composition of small local time intervals where the system is considered as monofractal.Now, the empirical data that will be used corresponds to 2 financial time series stated in Table 2.

Relationship of generalized Hurst exponent and the temporal Theil scaling exponent in empirical time series
The process to follow in this section consists of: 1. Take the closing price on each time series defined on a daily basis denoted as S t .
2. Carry out profiling of the time series filtering the days in which the empirical data have a non-zero return.
3. Estimate the logarithmic returns L t = ln(S t+1 ) − ln(S t ), and absolute log-returns, and volatilities of the log-returns using Eqs ( 33) and (34), respectively.Thus, the total time series is tripled.
4. Apply the diffusive algorithm for the original time series, absolute log-returns, and volatilities of the log-return.
5. Calculate the cumulative mean, and the Shannon index in each of the diffusive trajectory time series by adding new data for each time t.
6. Perform power law regression based on Eq (24) to find the temporal Theil scaling coefficient K TTS (t) and the temporal Theil scaling exponent α TTS (t).It is important to mention that in all the iterations of this process, one cost function is estimated on the adjustment of the temporal Theil scaling given by the coefficient of determination R 2 TTS .7. Estimate the generalized Hurst exponent using a multifractal method such as the MF-DFA method discussed in section 3.Then, to avoid the overestimate of the generalized Hurst exponent in short time series that it has already been shown in [65], a threshold value T MF is chosen such that H(q) is calculated if and only if t � T MF , where t represents the size of the cumulative time series at time t.
8. Filter the simulations where R 2 TTS exceeds a certain threshold value r TTS � 0 and t � T MF .Thus, we have that the selected data is from a time t � T MF which must be chosen appropriately since the diffusive trajectory algorithm is of the order Oðt 2 Þ. 9. Perform polynomial regression based on Eq (44) to find the coupling coefficients of the generalized Hurst exponent H(q) with the temporal Theil scaling exponent α TTS denoted by fb n ðqÞg W n¼1 .It is important to mention that in all the iterations of this process, two cost function is estimated on the adjustment of Eq (44) given by the coefficient of determination R 2 (q) and the p-th mean absolute error MAE p (q) defined in Eq (35).
10. Repeat the above process with different values of q, say W values, i.e. with q 2 fq k g W k¼1 .Finally, before going on to the results after having done this process, it is worth mentioning that the parameters selected were r TTS = 0.95, Regarding the estimation of the generalized Hurst exponent with the MF-DFA method, a library built in Python is used to make this calculation quite efficient [61].In addition, all the codes made are published in the Github [67].At last, to obtain the uncertainty in the parameters estimated by the least squares regression in some adjustments, the covariance matrix Cov-(�) is calculated with the vector of residuals of the fitẽ, that is, Covð�Þ / ẽẽ T .Then, the square root of the respective term on the diagonal of the covariance matrix is taken as the error in the fit parameters.
Figs 2 and 3 show the generalized Hurst exponent as a function of the temporal Theil scaling exponent for the given parameters.Hence, it is clear that in all cases the polynomial fit  works quite well to capture the fluctuations of the generalized Hurst exponent as a function of the temporal Theil scaling exponent.Specifically, the points represent the empirical data obtained on the 2 financial time series while the solid line represents the polynomial fit for the different values of q.Indeed, the smallest coefficient of determination R 2 (q) is 49.56% for the EURCOP = X with q = −2.Also, it is worth noting that the scales of all the figures respect that the generalized Hurst exponent is about 2 orders of magnitude larger than the temporal Theil scaling exponent.It is important to note that regardless of the case, the range covered by the temporal Theil scaling exponent is of the order of 10 −4 * 10 −3 .Furthermore, Table 3 shows the coefficients of determination R 2 (q) obtained in each case and also shows the highest value obtained from the regression coefficients, that is Hence, in Table 3 all the values B(q) are in the same order of magnitude around 10 8 and coincide with the coefficient of the polynomial with the highest degree, that is, a W TTS , which implies that the contribution of this term is really of the order of 10 8−4W * 10 8−3W , that is, 10 −8 * 10 −4 .
Finally, it is worth mentioning that in Figs 2 and 3, it is observed that the relationship between the generalized Hurst exponent and the temporal Theil scaling exponent is not a bijection between these two quantities due to the non-monotonic behavior of the curves, which becomes clearer when remembering the polynomial relationship of Eq (44).Nevertheless, it is important to remember that the LA-fBm is a local approximation, which implies that this non-monotonic behavior must be interpreted in terms of the time series having a multifractal behavior presenting moments of smaller or greater long-range correlation even when the temporal Theil scaling exponent increases.
Table 3. Parameters of the relationship between the generalized Hurst exponent H(q) and the temporal Theil scaling exponent α TTS obtained after the fit with the expression (36).Note that for currency ticker, = X is removed in its name.

Optimal selection of moments for the calculation of the generalized Hurst exponent
Now, noting the high orders of magnitude of these fit parameters in Table 3, it can be thought that b n (q) is associated with a measure of the most optimal value q with which to calculate the generalized Hurst exponent in a time series.To do this, from now on we set (for simplicity) W ¼ W þ 1, and observe that Eq (44) generates the following system of equations |ffl ffl ffl ffl ffl ffl ffl ffl fflffl {zffl ffl ffl ffl ffl ffl ffl ffl fflffl } |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } Hðq 1 ; q 2 ; . . .; q Wþ1 Þ ¼ Aðq 1 ; q 2 ; . . .; q Wþ1 Þã TTS : Therefore, it is worth remembering that the spectral norm of Aðq 1 ; . . .; q Wþ1 Þ evaluates the square root of the largest eigenvalue of the matrix A T A. This computation provides an approximation of the most significant eigenvector in the matrix, denoted as z ¼ ðzðq 1 Þ; zðq 2 Þ; . . .; zðq Wþ1 ÞÞ.Consequently, selecting the largest absolute component of z enables the estimation of q 2 q 1 ; . . .; q Wþ1 , which contributes the most to the computation of the generalized Hurst exponent given the temporal Theil scaling exponent α TTS (t).In simpler terms, the spectral norm of Aðq 1 ; . . .; q Wþ1 Þ and its corresponding eigenvector furnish an algorithm for identifying the optimal q-th moment of the probability distribution associated with an empirical time series to enhance the accuracy of generalized Hurst exponent estimation given the temporal Theil scaling exponent α TTS (t).In fact, it is important to remember that precisely one of the sources of multifractality in financial time series is the probability distribution associated with said time series rather than the temporal correlations between different time periods that can be destroyed by reordering randomly the time series as shown in [68].
It's important to emphasize that the q-th optimal moment algorithm is inapplicable to fractional Brownian motion.In the case of a monofractal system like fractional Brownian motion, selecting any arbitrary value for q would be sufficient to obtain the Hurst exponent H.In this context, constructing the matrix A from the expression (46) and considering Eq (31) results in a matrix of size 2 × 1, as only one q value and two coefficients are necessary for the linear regression (see Eq (31)).Consequently, the spectral norm would be computed for a 1 × 1 matrix, which is trivial for the q-th optimal moment selection algorithm.As a result, this approach would not reveal any significant difference in the system's behavior, unlike the two financial time series described in Table 2 as shown below.Now, the matrix A has been computed for each financial time series outlined in Table 2. Subsequently, Table 4 displays the spectral norm [69] alongside the associated eigenvector components.Indeed, the optimal value q 2 À 2; À 1; 1; 2 for ⌃DJI and EURCOP = X is determined as 1 and −2 for absolute log-returns and volatilities of log-returns, respectively.Specifically, in the case of EURCOP = X, components corresponding to q = −2 and q = 1 exhibit slight differences.This observation suggests that, for the time series EURCOP = X, minor fluctuations play a more crucial role in the system's behavior, whereas for ⌃DJI, larger fluctuations hold greater significance.This inference aligns with the known stability of EURCOP = X as a currency over the presented time period in Table 2, signifying its frequent occurrence of both positive and negative returns with relatively small values.Similarly, ⌃DJI, being a stock index that has experienced consistent growth during the study period in Table 2, reflects predominantly positive returns, with occasional minor negative returns compensated by substantial increases in the stock index's value.
Lastly, for the sake of comparing this method of selecting the q-th optimal moment value with other proposed methodologies, such as the MF-DFA method, it is noteworthy that in computing the generalized Hurst exponent using MF-DFA, it is customary to identify the q-th optimal moment where the overall detrended fluctuation functions Fðq; dÞ � d HðqÞ and Fðq þ 1; dÞ � d Hðqþ1Þ display significant similarity [31,61].Consequently, if the q-th moment is varied continuously, this convergence typically occurs (as a first approximation) whenever db n ðqÞ dq ¼ 0, for all n 2 f1; 2; . . .; Wg, since from Eq (44) it follows that for all η > 0. Hence, the optimal q-th moment value in the MF-DFA method is chosen as the one where the generalized Hurst exponent as a function of the temporal Theil scaling exponent does not vary much in its parameters fb n ðqÞg W n¼1 .Nonetheless, this approach may introduce bias into the matrix A by generating two rows with closely resembling values, making it challenging to discern the most optimal q-th moment value.Therefore, it is crucial to emphasize that the selection of the q-th optimal moment presented in this study involves identifying the most representative q-th moment value during the computation of the generalized Hurst exponent across a range of values q 1 ; q 2 ; . . .; q W .This selection process aims to capture the most significant fluctuations within a time series, whether positive or negative, as shown with the series of time of ⌃DJI and EURCOP = X.Table 4. Spectral norm and eigenvector associated to the matrix A used to estimate the optimal value q for each of the 2 financial time series defined in the Table 2 according to Eq (46).Note that for currency tickers, = X is removed in its name.

Conclusions
In summary, we establish a theoretical link between the Hurst exponent H and the temporal Theil scaling exponent α TTS through the multifractality partition function approach.Specifically, Section 4 introduces the generalized Shannon index SðqÞ, extending the Shannon index, while subsection 5 defines the partition function O TTS (q, δ) and the multifractal exponent of the generalized Shannon index β TTS (q) that relates with the generalized Hurst exponent H(q).For a fractional Brownian motion, this exponent β TTS (q) is expressed in Eq (31).Also, in Section 6, we validate the α TTS and H relationship for fractional Brownian motion.Multiple simulations with 512 time steps show a linear regression with positive λ TTS−MF and negative κ TTS−MF coupling constants.The R 2 is 84.96% for volatilities of log-returns, indicating a highquality fit.Finally, Section 7 applies the relationship between the temporal Theil scaling exponent α TTS and the generalized Hurst exponent H(q) to construct a selection algorithm for the optimal qth moment, obtaining q ¼ 1 and q ¼ À 2 values for ⌃DJI and EURCOP = X, respectively.The algorithm suggests that small fluctuations are more relevant for EURCOP = X, while large fluctuations dominate the behavior of ⌃DJI.All results are accessible on Github [67].
From a practical standpoint, it's crucial to note that while the findings in this article pertain to two specific financial time series, the algorithm for determining the optimal q-th moment may be applied to other financial datasets.This addresses two common challenges associated with employing the Hurst exponent in financial time series analysis.Firstly, it resolves the issue of selecting the most suitable q-th moment for estimating the generalized Hurst exponent.Secondly, it provides a solution for calculating the Hurst exponent in scenarios where data distribution is limited, i.e., when dealing with short time series.This is particularly significant as the relationship between α TTS and H(q) (see Eq (44)) implies that the Generalized Hurst exponent could be estimated by knowing the temporal Theil scaling exponent, not only in the case of a fractional Brownian motion.Also, with respect to this second problem, it is worth mentioning that there are other methods to estimate the Hurst exponent in short time series, for example, by maximum likelihood which has been verified in a fractional Brownian motion and an extension of this stochastic process known as fractional Brownian bridge [70].Nevertheless, since the comparison of the Hurst exponent in short time series using different methods is beyond the scope of this article, it is proposed as a future direction of work.

Fig 1 .
Fig 1. Temporal Theil scaling exponent as a function of the Hurst exponent for a fractional Brownian motion (fBm).(A) absolute value of fractional Brownian motion.(B) volatilities of the fractional Brownian motion.In all cases, the simulated data was obtained after averaging over N fBm = 200 simulations of T fBm = 512 time steps for each value of the Hurst exponent displayed.In addition, the error bars were taken with the standard deviations generated by the different simulations and the theoretical fitting to Eq (31) as a continuous line.https://doi.org/10.1371/journal.pone.0303252.g001

Fig 2 .
Fig 2. Generalized Hurst exponent H(q) as a function of temporal Theil scaling exponent α TTS (t) for the time series of the Dow Jones Industrial Average measured daily from January 3, 1992, to June 7, 2023.(A) absolute log-returns time series.(B) volatilities of the log-returns time series.In all cases, the empirical data is shown as points and the theoretical fitting to (44) with W = 4 as a line.https://doi.org/10.1371/journal.pone.0303252.g002

Fig 3 .
Fig 3. Generalized Hurst exponent H(q) as a function of temporal Theil scaling exponent α TTS (t) for the time series of the Euro to Colombian peso currency measured daily from January 2, 2003, to June 7, 2023.(A) absolute log-returns time series.(B) volatilities of the log-returns time series.In all cases, the empirical data is shown as points and the theoretical fitting to (44) with W = 4 as a line.https://doi.org/10.1371/journal.pone.0303252.g003

Table 2 . Stock market index and currency used to explore the relationship between temporal Theil scaling and generalized Hurst exponent from the generalized Shannon index approach.
Dates are placed in the ISO universal date format. https://doi.org/10.1371/journal.pone.0303252.t002